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Abstract 

Among the many methods available for the determination of inviscid fluxes across a sur- 
face of discontinuity, the flux-difference-splitting technique that employs Roe-averaged variables 
has been used extensively by the CFD community because of its simplicity and its ability to cap- 
ture shocks exactly. This method, originally developed for perfect gas flows, has since been 
extended to equilibrium as well as nonequilibrium flows. Determination of the Roe-averaged 
variables for the case of a perfect gas flow is a simple task; however, for thermal and chemical 
nonequilibrium flows, some of the variables are not uniquely defined. Methods available in the 
literature to determine these variables, seem to lack sound bases. The present paper describes a 
simple, yet accurate, method to determine all the variables for nonequilibrium flows in the Roe- 
average state. The basis for this method is the requirement that the Roe-averaged variables form 
a consistent set of thermodynamic variables. The present method satisfies the requirement that 
the square of the speed of sound be positive. 

Symbols 

A inviscid flux Jacobian 

a, b factors defined in Eqn. (5) 

c speed of sound 

Cpi specific heat at constant pressure for species i 

Cv,t,i specific heat at constant volume for translational and rotational energy for species i 

^ v; K i specific heat at constant volume for vibrational energy for species i 

E =e + (i? + v 2 )/2, total energy (sum of internal and kinetic energies) per unit mass 
e tr sum of translational, rotational, and formation energies per unit mass of the gas mixture 

e tr,i sum of translational, rotational, and formation energies per mole of species i 

e v vibrational energy per unit mass of the gas mixture 



e Vi vibrational energy per mole of species i 
F| G inviscid flux vectors in x- and y-directions 

fi constants, 1.5 and 2.5 for atomic and molecular particles, respectively 

H =E + p/p, total enthalpy of the gas mixture 
n x , n y x- and y-components of the unit outward normal to an interface 
p pressure 

q =(u 2 + v 2 )/2, kinetic energy 

R universal gas constant 

R, R' 1 left- and right-eigenvector matrices of A 

r 

T translational temperature 

Ty vibrational temperature 

T re f reference temperature, zero Kelvin 

t time 

U unknowns vector 

u, v velocity components in x- and y-directions, respectively 

W vector of production terms 

Wj molecular weight of species i 

Xj mass fraction of species i 

a i dp/dp. 

P dp/d(pE) 

Y specific heat ratio in ideal gas 

A hfi heat of formation of species i at reference temperature 

A vector of eigenvalues of matrix A 

<t> dp/d[ pe y \ 

p density of the gas mixture 
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p ( - partial density of species i 

Subscripts: 
e electrons 

i, j indices for chemical species 

n normal to the interface 

L left side of the interface 

R right side of the interface 

t translational mode 

V vibrational mode 


Introduction 

In the numerical simulation of high speed flows, the flux-difference-splitting techniques 
have been successfully employed in the determination of the inviscid fluxes across a surface of 
discontinuity. Among these techniques, the one that uses the so-called Roe-averaged variables is 
popular among the researchers because of its simplicity and ability to capture shocks exactly. 

This is true for the case of perfect gases for which Roe [1] derived expressions for the necessary 
variables in the Roe-average state. The procedure has since been extended to equilibrium flows 
(see [2]-[5]), and to chemical and thermal nonequilibrium flows (see [6]-[8]). However, in all 
these extension the simplicity associated with a perfect gas case is lost. In addition, some of the 
Roe-averaged variables are not defined uniquely. This led to several methods for the determina- 
tion of these variables, some of which are based on computational experience. However, these 
methods lacked sound bases in their approach. 

The present paper addresses this issue, and describes a method for the determination of the 
Roe-averaged variables for thermal and chemical nonequilibrium flows. The method is simple, 
accurate, and requires no approximations. The basis for the present method is a requirement that 
the variables in the Roe-average state form a consistent set of thermodynamic variables. For the 
case of an air model with charged particles, the procedure requires the determination of a vibra- 
tional temperature for given values of vibrational energy and mixture composition. This, how- 
ever, is a routine computation. For the case of a gas mixture with only inert particles, the method 
gets considerably simplified. Finally, it is important that any such method yield physically mean- 
ingful results in all cases, i.e., the square of speed of sound in the Roe-average state must be posi- 
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tive. The present method satisfies this requirement for all cases. 


Analysis 

Consider the Euler equations for two-dimensional flows written in conservation form: 


dU 3F 3G 
dt dx dy 


= 0 


( 1 ) 


where 


U = [p, pu, pv, p£] T 

F = [p u,p + pu 2 , p MV, P uH] T (2) 

G = [pv, pvu,p + pv 2 , pvH] 


We consider two-dimensional flows only for simplicity; the present analysis carries to three- 
dimensional flows without any loss of generality. If a computational domain is discretized into 
cells, and the unknowns are assumed to be discontinuous across the interface between two adja- 
cent cells, then a question arises as to what is the flux across these surfaces of discontinuity? This 
is a classical Riemann problem for which several approximate solvers are available. Among 
these solvers, the one that uses Roe-averaged variables has been widely used by researchers in the 
CFD community. Determination of the Roe-averaged variables for the perfect gas case is simple, 
and these variables satisfy the jump conditions across normal shocks exactly. Adopting this 
approach, the numerical flux across a surface of discontinuity is written as 

F » - [^+r„« + lA|(O i ,-U l )]/2 (3) 


where and U L are the unknowns in the left and the right cells, respectively, and ¥ nL and F n/? 
are the fluxes in the left and right cells normal to the common surface, respectively. The matrix 
A is the Jacobian of the inviscid fluxes, and is evaluated at the Roe-averaged state U , which in 
turn is a function of U L and . The matrix |a| = r|a|r where the matrices R and R 1 are 
the left and the right eigenvector matrices of the matrix A . Roe [1] identified the properties this 
matrix A should satisfy as the following: 

(i) The matrix A constitutes a linear mapping from U to F . 

(ii) As U L — » Ufl — » U, A (U L , U^) — > A (U) , where A = 3F n /3U . 

(Hi) For any U L and , A (U L , U^) • (U L - U^) = F L -F R . 
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(jv) The eigenvalues of A are linearly independent. 

Roe determined the unknowns U , and hence the matrix A that satisfies these properties. The 
unknowns that determine the matrix A , its eigenvalues, and its eigenvectors are given by 

u = au L + bu R v = av L + bv R H = aH L + bH R (4) 

where 

* a = JPr] = 17 ( 1 +0 b = J^R / [J^L + JPR~] = r/ (' +r ) ( 5 ) 

Note that r = Jp R / and a + b = 1 . The speed of sound c is given by 

c 2 = (y-l) (H-q) (6) 

where q = (u 2 + v 2 )/ 2. Computing the variables u, v , and H using Eqn. (4) has since become 
known as Roe-averaging. Flux difference splitting technique using Roe-averaged values has been 
successfully used by many researchers in solving high speed flows. 

Nonequilibrium Flows 

Consider the governing equations for thermal and chemical nonequilibrium flows in two- 
dimensions. The number of governing equations isN = n s + 2 + n e , where n s is the number of 
chemical species considered in the chemistry model, and n e is the number of energy equations. 
The other two equations come from conservation of *- and y-momentum. Gnoffo [6 & 9] consid- 
ers a 11 -species air model with two temperatures - translational and vibrational. In determining 
the flux across a surface of discontinuity, he employs the flux-difference-splitting technique. 
However, explicit expression for the variables in the Roe-average state are available for only 
some of the variables. In his analysis, Gnoffo [6] uses these expressions, and makes approxima- 
tions for the remaining variables, based on his computational experience. Other approximate 
methods for the determination of these variables are available, (see Liu & Vinokur, [7], for exam- 
ple). In the present paper we examine Gnoffo’s multi-species, two-temperature model, and 
present a method for the determination of the Roe-averages for all of the necessary variables with- 
out making any approximations. 

As in the ideal gas case, we consider two-dimensional flows for simplicity. The governing 
equations for nonequilibrium flow of a multi-species, two-temperature gas written in conservation 
form as described by Gnoffo [9] are as follows: 
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(7) 


8U 3F 3G 

dt dx dy 


= W 


where 


U = [p p p u, pv, p E, pe v ] T 
r 2 -i t 

F = [p^, p + pu , puv, puH, p«e v J (8) 

G = [pj.v, pvu,p + pv 2 , puH, pveJ T 

Only the convective terms are considered. The quantity W Eqn. (7) is a vector of the production 
terms, and does not enter into the present analysis. Most of the analysis described here follows 
Gnoffo’s analysis with minor changes in the notations. 

There are n s mass conservation equations - one for each chemical species considered, two 
momentum conservation equations, and two energy conservation equations - one for total energy 
and one for vibrational energy. The term H is the total enthalpy, and is equal to E + p/p where E 
is the total energy given by e (r + e v + q . The quantity e tr is the sum of translational, rotational, 
and formation energies of molecular and atomic particles, e y is the sum of vibrational, electronic, 
and electron translational energies, and q = (u 2 + v 2 ) /2 is the kinetic energy. In the two-tem- 
perature model the translational temperature, T , describes e tr , and vibrational temperature, T v , 
describes e y . Henceforth, we shall use the terms ‘vibrational energy’ to represent the sum of 
vibrational, electronic, and electron translational energies. Note that the temperatures do not 
appear explicitly in the governing equations. Also note that in the vibrational energy equation, 
the term corresponding to the work done by the electron pressure gradient has been moved over to 
the right hand side of the governing equations and lumped with the production terms. This has 
lead to some simplifications in the analysis. 

Before proceeding to the determination of the Roe-averaged variables and computing the 
fluxes, it is useful to review the procedure for the determination of the energies. In the two-tem- 
perature model, the translational and rotational modes are assumed to be fully excited, and hence 
the heat capacities of these modes are independent of the temperature, and the internal energy is a 
linear function of the translational temperature, i.e., 
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(9) 


T 

<v.<m = J c v ,,,^r)dr + Ah fi 

T 

1 ref 

Without loss of generality we assume the reference temperature, T ,, to be equal to zero Kelvin. 
Note that C v r . are constants (independent of temperature), and are equal to 2.5 R/W i for molec- 
ular species and 1.5 R/W { for atomic species. The quantity Ah^ { is the enthalpy of formation of 
species i at the reference temperature. 

Unlike the translational-rotational modes, the vibrational and electronic energy modes are 
assumed to be only partially excited; hence the vibrational heat capacity depends on the vibra- 
tional temperature. The heat capacity for the vibrational mode of species i can be evaluated by 
utilizing the curve fits for the total heat capacity evaluated at T y , and subtracting out the constant 
contribution from translational and rotational modes as follows; 

C v , v, i ( = C v , i ( T v ) - c v t i ( 1 0) 

Note that C v • (T v ) = C p y ( T y ) - R/W r Curve fits for C p ( . in terms of T may be found in the 
literature (see e.g., Prabhu and Erickson [10]). Integrating C v v . from T ref (equal to zero 
Kelvin) to the required temperature T y we obtain e y ., i.e., 

T v 

e v ,i(T v ) = J C v V I (T')dr ( 11 ) 

T „< 

The energies per unit mass of the mixture are determined as follows: 

n s - 1 ", 

e tr ~ X e ir, i X i e V ~ X e r,i X i (12) 

i = 1 i = 1 

Note that the summation for translational energy is over all the species except the electrons; the 
electron translational energy is included in the vibrational energy. Equations (9) - (12) are rela- 
tions between the two energies, their respective temperatures, and the mixture composition. For 
given values of the energies and the mixture composition, the translational and vibrational tem- 
peratures can be determined. Similarly, if the temperatures and the mixture composition are 
known, the energies can be computed. These are routine computations. 

The component of flux normal to the common side between two cells is denoted by F . 

J n 

The Jacobian matrix A = 9F n /8U can be evaluated (see Gnoffo [9]), and is given by 
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( 13 ) 


U (8..-JC.) 

xn 

x.n 

0 

0 

IJ V 

1 X 

* y 

°A ~ Uu 

(1-P) un x + U 

( 1 - P) vn x - V 

K 


a n -Uv 
j y 

(1-P) un y + V 

(1-p) un y +U 

P", 

<1 >«y 

(a j-H)U 

Hn x -$Uu 

Hn x - pt/v 

d+P)t/ 

<\>u 

-Ue v 

e V*x 

e V n y 

0 

u 


The eigenvalues of A can be determined, and are given by 

A = (U repeated (n,+ l) times, U + c,U-c, U ) (14) 

Subscript i refers to species i and row i, and subscript j refers to species j and column j, and 8^ is 
the Kronecker delta. The quantities U = un x + vn y and V = vn x — un y are the velocity compo- 
nents normal and parallel to the cell interface respectively, n x and n y are the *- and y-components 
of the unit outward pointing normal to the interface respectively, * ( . are the mass fractions of the 
constituent species, and c is the speed of sound defined by 

", 

c 2 = I a,*, + P (H- 2q) + <|>e v (15) 

I = 1 

The quantities a { ,i = 1, n s , P, and <f> are the partial derivatives of pressure, defined as fol- 
lows: 


dp RT Q a 

~ dp~. - ~W~ + ~ $ e i ~ i 


n s - 1 

6 = dp _ ^ v Jj. 
dp E~ C v , £ W t . 


(16) 


(17) 


4, = *- 

dpey C v y 



(18) 


In Eqn. (16), T = T for heavy particles (atoms, molecules, and ions), and T 
trons. Note that the quantities C v t and C v v are given by 

C v,t - X X i C v, t , i C v, V = X X i C v, V, i 

i=l i=l 


T v for the elec- 


(19) 


As noted earlier, C v , . = 2.5 R/W i for molecular particles and 1.5 R/W i for atomic particles 
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making C v t dependant only on the mixture composition. Consequently, P also depends only on 
the mixture composition, and is easily determined if the mixture composition is known. The 
quantity <J) , however, depends on the mixture composition and C v v , which is a function of the 
vibrational temperature. 

The left- and right-eigenvector matrices of A corresponding to the eigenvalues in Eqn. 
(14) are denoted by R and R ' , and are given by 


r 2 

c o . . - a .x . 
V J ' 

-V 


R = 


a.- Uc 
J 

a . + Uc 
J 

-a.e 
J V 


5 .. 

V 


c 2 R 1 = 


[ 2 - <r] /p 


pMJt. 

1 

pvx. 

i 

~V x i 



-rt 

n 

0 

0 


y 

cn - B u 

y K 

X 

cn - Bv 

y 

P 

<t> 


-cn - B« 

-cn - Bv 

> ? 

P 

<P 


ftue y 

$ ve V 

-P<V 

2 - §e y 


0 

X . 
1 

x. 

1 

0 

2 

~C Hy 

[u + cnj /2 

[u-cn^ /2 

0 

C n x 

[ v + cn 1/2 

[v-cn ] /2 

0 

c 2 V 

y 

[ H + cU ] /2 

y 

[H-cU] /2 

-VP 

0 

e/2 

v /2 

1 


( 20 ) 


( 21 ) 


The matrices R, R , and the eigenvalues of A defined by the Roe-averaged variables 
(denoted by an overline) are required for the determination of the flux. In order to determine the 
Roe-averaged variables, we start with the following necessary condition that A must satisfy: 


A(F„) = AA (U) ( 22 ) 

where A( ) = ( ) R - ( ) L . We are seeking a set of variables that satisfies the above equation 
exactly. This equation has N components corresponding to the N governing equations. Of these 
N components, the first n s components (corresponding to the species mass conservation equa- 
tions) are satisfied, if we define 


x i = a (x { ) + b (x^) u = au L + bu R v = av L + bv R (23) 

where, as in the ideal gas case, the quantities a and b are functions of densities in the left and 
right cells, and are defined in Eqn. (5). Note that 
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(24) 


X * i = °X (*«■)* = <*+*> = 1 


i = 1 


I = 1 


/= 1 


This is a necessary condition if x it i - 1 , . . . , 5 are to retain the property of the mass fraction. 
Next, the components of Eqn. (22) corresponding to the momentum equations are satisfied pro- 
vided the condition 


A(p) = X (« J -P?)A(p i ) + pA(pe) + <DA(peJ 


(25) 


I = 1 


is satisfied. Note that e = e tr + e y . Finally, the components of Eqn. (22) corresponding to the 
total energy and the vibrational energy equations are satisfied, if we define 


H = aH L + bH R 


e v ~ ® ( e v) , + h (^v) 


(26) 


V' u R 

Thus, we have determined the unknowns x-, u, v, H , and e y in terms of the corresponding quanti- 
ties in the left and right cells, and the density ratio across the interface. Of these, the relations for 
u, v , and H are the same as in the ideal gas case. A condition similar to Eqn. (25) for the ideal 
gas case is satisfied identically. 

It remains to determine the quantities a„ p , and <j> which are not defined explicitly. The 
only condition they must satisfy is the Eqn. (25). We describe here a simple method to compute 
these quantities. The basis for the present method is the requirement that these unknowns 
together with x i and e y form a consistent set of thermodynamic variables. We consider two 
cases - (1) a gas mixture containing charged particles and electrons, and (2) a gas mixture of only 
neutral particles. The first case is treated in detail. The second case is much simpler, and the 
results follow from the those of the first case. 

Case 1: Gas mixture of charged particles and electrons:- 

As noted earlier, p depends only on the mixture composition. Hence, we require that the 
Roe-averaged quantity, P , depend only on the mixture composition at the Roe-average state, and 
be determined using a relation similar to Eqn. (19). Accordingly, we write 
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where, as in Eqn. (21), we note 


n, - 1 


n.-l 


C v-' ~ Z = Z (28) 

i = 1 i = 1 

Hence, p can be readily computed. Following a similar reasoning, we write for <j> as in Eqn. (20) 

(29) 


$ = 


r W 

c v, V e 


Again, as in Eqn. (21) we note 


Cv, v - X v, j (Tv) (30) 

i = l 

The value of C v v ,• , however, must be evaluated at an appropriate vibrational temperature, T v . 
Since e v and the mixture composition x ( are known, we require that the appropriate vibrational 
temperature, T v , be consistent with e v and x i . This temperature can be determined following a 
routine computation. Once Ty is determined, (|> can be evaluated. 

Next, we proceed to the determination of a, . Following the argument of the previous 

paragraph, we require a, to have the same functional relation with the other variables as in Eqn. 
(18), and write 


“ RT — - 

a * = (3i) 

In this equation, all the quantities except e ( and T are known. For electrons, T = T v , which is 

known, and e. = e V i = 1.5 RT V /W r Hence a, for electrons is readily obtained, and may be 
expressed as follows: 

a e = KT v [l-1.5(P + V]/W e + fe 
For the heavy particles we have T = T\ hence we write 


( 32 ) 



~ e i = *tr, i + ~ e v,i = fi RT/ w i + • + e v ■ (33) 

This equation together with Eqn. (31) provides a linear relation between a, and T , which may be 
written as follows: 

a i = *^[1 ~f$] /w i + P (9 - A/iy- ,.) - (p + <j>) f (34) 

This equation alone is not sufficient to solve for T (or a,- ). Moreover, we need to satisfy the con- 
dition in Eqn. (25). Hence we combine Eqns. (25), (32), and (34), and obtain the following: 

n,-\ 

A(p) = R?'£ [(l-/ I P)/W r PA/, /ir (P + ^)e v> jA(p.) 

« = i (35) 

+ (“ e ~P?)A(p e ) + pA(pe) +^A(pe y ) 

Solving for the only unknown T in this equation, we find 

A(p) - PA(pe) -<|>A(p<r v ) -<x c A(p e ) 

< n - 1 

5 r 

+ ^ [pAA /f + (P + ^)« y .]A(p.)-/tr v [l-1.5(P + ^)]A(p e )/W e 


R £ A(p.) [1 
i= 1 

Once T is determined, a, for the heavy particles can be computed using Eqn. (34). Although T 
does not occur explicitly in any of the matrices used in the determination of the fluxes, it is desir- 
able that the value of T computed by Eqn. (36) is positive. In order to demonstrate this, we sim- 
plify Eqn. (36) and write it in the following form: 

f = aT L + bT R~ 


(MljMpV- X^Atp.) +« A[p^r r fyJ]/w (37) 

irn 

s 

R X [ 1 -/ I P]A(P I )/W'. 
i = 1 

It has not been possible to further simplify Eqn. (37) so as to demonstrate that the values of f 
given by this equation are always positive. However, on the right-hand-side of this equation, the 
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first two terms are positive and large, and the quantities in the numerator of the third term, being 
differences of similar quantities, are small. In all the numerical exercises we have performed, the 
third term has always been only a small fraction of the first two, and the value of T has been a 
positive quantity. In the freestream this term is indeterminate, and is set equal to zero. Once T is 
determined, a, for the heavy particles can be computed using Eqn. (34). 

The Roe-averaged speed of sound may be determined using the following: 

", 

c 2 = £ <*,*,.+ p [H-2q] +<t>e v (38) 

i = i 

Thus, all the Roe-averaged values of the quantities ex,, i = 1, n s , p, <j) , and c required to com- 
pute the matrices A, R , and R are determined. These quantities satisfy the necessary condition 
given by Eqn. (25) exactly. Together with x., e y , T, and T v , they form a consistent set of thermo- 
dynamic variables, which is the basis for the method of determining these quantities described 
here. The method is simple, and makes no approximations. The only additional computation 
required is the determination of the vibrational temperature, T v for the given values of vibrational 
energy, e v , and the mixture composition, x. ; however, this is a routine computation. 

It is necessary that the value of the square of the speed of sound computed by Eqn. (38) be 
. . 2 

positive. To demonstrate that c computed by Eqn. (38) is always positive, we start with the def- 
inition of H 


H = aH L + bH R = ah L + bh R + q + 8 (39) 

where 8 = ~^[_(u R - u L ) 2 + (v R - v L ) 2 ] > 0 . Next, we substitute h = e {r + e v + p/p in the 
above equation, and obtain the following: 

H = e^ r + e y + p/p + q + 8 (40) 

where e tr = a (e tr ) L + b (e tr ) R and p/p = a (p/p) L + b (p/p) R . 

Next, in Eqn. (38), we substitute for a, from Eqns. (32) and (34), and for H from Eqn. 
(40), and simplify to obtain 

c 2 = P[£+ (p7p) +8] +RT v x/W e (41) 

where e = ae L + be R , and 
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( 42 ) 


71,-1 

e = RT £ f i x/W i 
i = 1 

All the terms on the right hand side of equation (41) are greater than or equal to zero, proving that 
the square of the speed of sound, c , is always positive. 

Case 2: Gas mixture of only neutral particles 

The procedure described above is simplified when electrons are absent in the gas mixture. 
The governing equations are the same as in Eqns. (7) and (8), and the thermodynamic relations 
given in Eqns. (9), (10), (11), (12), and (14)-(19) as well as the matrices A, R , and R 1 given in 
Eqns. (13), (20), and (21) are also applicable. 

The Roe-averaged quantities defined in Eqn. (23), (24), and (26) are applicable without 
any changes. The necessary condition given in Eqn. (25) remains the same as well. Equation 
(27) for the quantity p is applicable. However, since electrons are absent, the quantity <j> is given 
by 


<t> = P (43) 

and there is no need to determine the quantity Ty • With P + <j> — 0 , the expression for cq , given 
in Eqn. (34), can be written as follows: 

a i = RT[ 1 -/;.p] / W ( . + p (q - A h f ,.) (44) 

The temperature, T, is determined following a similar procedure, and is given by 


T = 


A ( p ) - pA ( pe (r ) + p j A (p.) A h f . 


i= 1 


n -1 


[l-/ / P]A(p.)/W. 


(45) 


i = 1 


This equation for T can be simplified and may be written in the following elegant form: 


T=aT L + bT R (46) 

This also demonstrates that the temperature T is always positive, and lies in between T. and T 
After computing T, a, can be determined using Eqn. (44). 

The speed of sound given in Eqn. (38) is still applicable and can be written as 
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( 47 ) 


n , 

^ = '5', &i x i + P [7/ — 2q — ey] 

« = l 

This can be further simplified, and written in the following form: 

c 2 = p (e + (pTp) + 8) (48 ) 

again proving that the square of speed of sound determined by this method is always positive. 

Conclusion 

A new method is presented to determine the Roe-averaged variables required to determine 
the inviscid fluxes in chemical and thermal nonequilibrium flows. These variables form a consis- 
tent set of thermodynamic variables, and satisfy the necessary condition. The method is simple, 
and does not make any approximations. The only additional computation required in determining 
the Roe-averaged variables for gas mixtures with electrons, is the computation of vibrational tem- 
perature for given values of the vibrational energy and the mixture composition. This, however, 

is only a routine computation. A logical requirement that the square of the speed of sound be pos- 
itive is satisfied. 
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